##########################################
#Author:Naresh Prodduturi
#Code:post Process for genotyping hazard rank boot strap methos
##########################################
#readin the calculated AUC ranks for all the iterations
args <- commandArgs(TRUE)
library(survival)
AUCinput<-read.table("top_genoexpression_AUC.txt",sep="\t",head=T,quote = "\"")
rows<-rownames(AUCinput)

#readind the original genotyping file
genotype<-read.table("top_input_expr",sep="\t",head=T)
number_top<-20
finalhazardrank<-matrix(nrow=number_top,ncol=ncol(genotype))
colnames(finalhazardrank)<-colnames(genotype)
maintop20<-c()

#loop over the boot strap iteration i1
i1<-1
#readin the martingale order
ordersamples<-read.table("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/ordered_sampes_Martingale.txt")
ordersamples<-as.vector(ordersamples[,1])
ordersamples<-gsub("-",".",ordersamples)
s1<-seq(1:207)
s2<-sample(s1, size=length(s1), replace = FALSE, prob = NULL)
ordersamples<-ordersamples[s2]
# reading the clinical file
	clin<-read.delim("/data1/bsi/BORA_processing/devel/eqtl/parallelize_genotypingsurvival/clinical" ,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="",header=TRUE,col.names=c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE"),as.is=c(TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),colClasses=c("character","numeric","factor","numeric","numeric","numeric","integer","integer","integer","factor"))
	clin$SUBJECT<-gsub("-",".",clin$SUBJECT)
	samplesize<-20
	#num<-as.integer((length(geno[1,]))/samplesize)
	
	#sample size 20
	num<-20
	
while(i1<ncol(AUCinput)+1)
{
	#creating uniqure names to store for hazard rank ,AUC and top20 snplist each iteration 
	hazardname<-paste("hazardname",i1,sep="")
	aucname<-paste("auc",i1,sep="")
	snplist<-paste("snp",i1,sep="")
	
	#sorting AUC rank for that column and storing index in list vector
	list<-rev(sort.list(AUCinput[,i1]))
	
	# selecting top 20
	top20<-rows[list[1:number_top]]
	maintop20<-c(maintop20,top20)
		
	#code start for hazard rank and AUC
	#preparing top 20 snps for hazard rank calculations and AUC 
	geno<-genotype[list[1:number_top],]
	i<-1
	
	#looping through each snp and calcuating hazard rank and AUC
	while(i<nrow(geno)+1)
	{
		#excluding the annotation
		geno1<-geno[i,]

		#code for craeting the iteration blocks sub and main main part_1_19 sub 20
		d<-i1
		s<-d;
		sub<-c()
		main<-c()
		if(d > 1)
		{
			main<-seq(1,d-1,1)
		}
		nu<-1
		while(s<length(geno1)+1)
		{
			sub[nu]<-s
			sequ<-c()		
			if(s+1 <= length(geno1))
			{	
				sequ<-seq(s+1,s+num-1,1)
				if(s+num-1 > length(geno1))
				{
					sequ<-seq(s+1,length(geno1),1)
				}
			}
			main<-c(main,sequ)
			nu<-nu+1	
			s<-s+num
		}
		#code for creating iteration block ends
		genomain<-geno1[main]
		genosub<-geno1[sub]
		dt <- merge(clin,t(genomain),by.y="row.names",by.x="SUBJECT")
		colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		
		j=1
		order_samp<-c()
		while(j<length(ordersamples))
		{
			order_samp<-c(order_samp,which(dt[,1]==ordersamples[j]))
			j<-j+1
		}
		order_samp<-c(order_samp,which(dt[,1]==ordersamples[length(ordersamples)]))
		dt<-dt[order_samp,]
		
	
		#create survival model
		fit1 <- coxph(Surv(OS,ALIVEOS) ~ EXPR+AGE,data=dt)
		dt1 <-merge(clin,t(geno1),by.y="row.names",by.x="SUBJECT")
		colnames(dt1) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
		
		j=1
		order_samp<-c()
		while(j<length(ordersamples))
		{
			order_samp<-c(order_samp,which(dt1[,1]==ordersamples[j]))
			j<-j+1
		}
		order_samp<-c(order_samp,which(dt[,1]==ordersamples[length(ordersamples)]))
		dt1<-dt1[order_samp,]
		
		pred <- predict(fit1, dt1,type="lp")
		#code for sorting start
			order_rank<-order(pred)
                	sorted_ranks<-seq(1,length(order_rank),1)
                	order_of_index<-order(order_rank)
                	wanted_ranks<-sorted_ranks[order_of_index]
		#code for sorting end
		finalhazardrank[i,sub]<- wanted_ranks[sub]
		i<-i+1
	}
	i1<-i1+1
}
i<-1
AUCrank<-c()
while(i<nrow(finalhazardrank)+1)
{
	dt <- merge(clin,finalhazardrank[i,],by.y="row.names",by.x="SUBJECT")
	colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","RANK")
	surv<-survConcordance(Surv(OS,ALIVEOS)~RANK,data=dt)
	AUC<-surv$concordance
	AUCrank<-c(AUCrank,surv$concordance)
	i<-i+1
}
#write.table(AUCrank,"finalAUCrank.txt",quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
write.table(AUCrank,args[1],quote=FALSE,sep="\t",col.names=TRUE,row.names=TRUE)
